 setwd("E:/1.甲基化分析/ASM/ASM_snp-onlyWGS")
file=read.table("E:/2_5hmc_yjp_bam/5hmc_file/2_5hmc_yjp_bam/DMR_edgeR_noRMS/6v6p_8v8p/chipseq_encode/motif_overlap_TF_chipseq_meta.txt",head=T,sep="\t")
file1=file[file$BF_in_DC>10&file$metap_in_DC_CC<0.05,]


ID=as.character(file1$unitID)
sel=list.files(pattern="snp")
ID=c("chr18:5289018:5297052")
pick=function(df,chrs,start,end){
df=df[df$chrom==chrs,]
dfs=df[df$position>start&df$position<end,]
return(dfs)
}


for (i in 2:length(sel)){
tmp=read.table(sel[i],head=T,sep="\t")
tmp$unitID=paste(tmp$chrom,tmp$position,sep=":")
tmp1=tmp[1,]
#tmp1$tartget="none"
for(j in ID){
chrs=unlist(strsplit(j,":"))[1]
start=unlist(strsplit(j,":"))[2]
end=unlist(strsplit(j,":"))[3]
tmps=pick(tmp,chrs,start,end)
#tmps$tartget=j
tmp1=rbind(tmp1,tmps)
}
tmp1=tmp1[-1,]
tmp1$tumor_var_freq=as.numeric(gsub("%","",tmp1$tumor_var_freq))/100
tmp1$normal_var_freq=as.numeric(gsub("%","",tmp1$normal_var_freq))/100
tmp2=data.frame(tmp1[,5:7],tmp1[,9:11],exitS=tmp1$exits,unitID=tmp1$unitID)
fn=unlist(strsplit(gsub(".snp","",sel[i]),"_"))
names(tmp2)=c(paste0(fn[1],c("_reads1","_reads2","_var_freq")),paste0(fn[2],c("_reads1","_reads2","_var_freq")),paste0(fn[1],"_",fn[2],"_exits"),"unitID")
final=merge(final,tmp2,by="unitID",all=T)
}

###find some loci
ID=c("chr18:5297603")
for (i in 2:length(sel)){
tmp=read.table(sel[i],head=T,sep="\t")
tmp$unitID=paste(tmp$chrom,tmp$position,sep=":")
tmp1=tmp[tmp$unitID==ID,]
tmp1$tumor_var_freq=as.numeric(gsub("%","",tmp1$tumor_var_freq))/100
tmp1$normal_var_freq=as.numeric(gsub("%","",tmp1$normal_var_freq))/100
tmp2=data.frame(tmp1[,5:7],tmp1[,9:11],exitS=tmp1$exits,unitID=tmp1$unitID)
fn=unlist(strsplit(gsub(".snp","",sel[i]),"_"))
names(tmp2)=c(paste0(fn[1],c("_reads1","_reads2","_var_freq")),paste0(fn[2],c("_reads1","_reads2","_var_freq")),paste0(fn[1],"_",fn[2],"_exits"),"unitID")
final=merge(final,tmp2,by="unitID",all=T)
}